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1  Executive  summary 

In  the  light  of 

-  the  seriousness  of  the  noise- induced  hearing  damage, 

-  the  lack  of  fundamental  knowledge  regarding  transmission  of  acoustic 
energy  through  non-airborne  pathways  to  the  cochlea,  and, 

-  difficulties  associated  with  designing  experiments  reproducing  realistic 
conditions  of  sound  propagation  through  the  human  head, 

our  objective  has  been  to  develop  rigorous  algorithms  and  to  construct  corre¬ 
sponding  efficient  software  tools  allowing  high  fidelity  numerical  simulations 
of  sound  wave  propagation  through  biological  tissues- 

As  the  main  achievements  of  our  Phase  I  work  we  consider: 

■  development  of  a  set  of  algorithms  (including  non- lossy,  error- 
controlled  matrix  compression  techniques)  for  general,  Le.,  variable- 
density,  volumetric  equations  of  acoustics;  and 

«  development  of  a  solution  technique  overcoming  ill-conditioning  diffi¬ 
culties  arising  in  large  density  contrast  problems. 

We  believe  that,  in  this  way,  we  contributed  a  practically  useful  and  efficient 
approach  for  solving  a  largo  class  of  problems  of  acoustic  wave  propagation 
through  inhomogeneous  media,  including  propagation  of  sound  through  the 
human  body. 

To  execute  the  proposed  effort,  we  combined  capabilities  and  expertise  in 
the  areas  of  ; 

-  development  of  fast  solvers  for  realistic,  parallel,  large  scale  simulations 
(Monopolc  Research), 

-  generation  of  accurate  geometry  representations  (UT  Austin)  and, 

-  identification  and  execution  of  relevant  acoustic  wave  experiments  (UT 
Austin). 


1 


2  Objective  and  brief  approach  description 


To  reiterate,  the  objective  of  our  effort  was 

to  develop  powerful  software  tools  and  to  perform  high  fidelity 
simulations  which  would  allow  identification  and  understanding 
of  relevant  bioacoustic  and  psychoacoustic  mechanisms  responsi¬ 
ble  for  the  transmission  of  acoustic  energy  through  non- airborne 
pathways  to  the  cochlea . 

As  our  tool  for  simulating  sound  propagation  in  the  human  body,  we 
chose  an  approach  based  on  volumetric  integral  equations.  The  well  known 
advantages  of  the  integral-equation  formulation  include  its  high  numerical 
accuracy,  stability,  and  exact  treatment  of  the  boundary  conditions  at  infin¬ 
ity  (and,  hence,  no  need  of  discretizing  any  regions  of  the  space  surrounding 
the  scattercr). 

The  main  reason  why  integral  equations  had  not  been  more  commonly 
used  in  acoustics  is  their  high  computational  cost:  a  straightforward  integral- 
equation  discretization  Leads  to  the  costs  of  at  least  G(N2)  for  iterative  and 
O(N^)  for  the  direct  solution  methods.  However,  adaptation  of  presently 
available  matrix  compression  and  fast  solution  techniques  [1,  2]  to  acoustics, 
as  described  in  Section  3.3,  reduces  that  cost  to  approximately  0(N },  and 
thus  renders  large  integral-equation  problems  tractable. 

Another  issue  arising  in  applying  the  integral-equation  formulation  to 
sound  scattering  on  a  body  surrounded  by  air  is  the  necessity  of  taking  into 
account  the  spatial  variability  of  the  medium  density ,  as  the  ratio  of  a 
biological  tissue  density  to  that  of  air  is  about  1000.  Our  development  of  a 
two-stage  solution  technique  constitutes,  to  our  knowledge,  the  first  practical 
numerical  implementation  which  overcomes  ill-conditioning  caused  by  large 
density  contrasts,  and  results  in  a  well-posed  problem.  1 

In  what  follows  we  briefly  describe  the  main  results  and  accomplishments 
of  our  Phase  I  work. 

3  Summary  of  Phase  I  results 

3.1  Formulation  and  implementation  of  an  integral-equation 
approach  in  acoustics 

We  present  here  the  motivation  of  our  choice  of  the  integral  equations  as 
the  mathematical  modeling  tool,  and  give  a  short  description  of  our  integral 

1  We  note  that  much  of  the  existing  work  on  ^acoustics  integral  equations”  has  been 
concerned  exclusively  with  constant- density  problems. 
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equations  and  their  discretization.  Further  technical  details  are  given  in 
Appendices  A  and  B. 


Volumetric  integral  equations.  Volumetric  integral  equations  for  the 
pressure  field  -  known  generally  as  the  Lippmann-Schwinger  (L-S)  equations 
-  can  be  derived  from  the  differential  equation  of  acoustics  split  into  two 
terms:  one  describing  wave  propagation  in  the  background  medium  (in  our 
case,  air)  and  the  other  wave  interaction  with  the  scatterer  (here,  the  human 
head),  occupying  the  spatial  region  ft. 

The  resulting  L-S  equation  [3]  for  the  pressure  field  p(r)  is 


P<mc){ r)  =  p(r)  +  J  dV  | k? g(r  —  r')  (l  -  P(r') 

n  ^ 

-  [vr,  ff(r  -  r')]  ■  (l  -  vr,p(r')|, 


(3.1) 


where  p(inQ)  is  the  incident  wave,  satisfying  the  Helmholtz  equation  in  the 
background  medium,  and  g  is  the  corresponding  Green  function, 


0(r) 


47r|r| 


(3-2) 


The  details  of  the  derivation  arc  given  in  Appendix  A. 


Discretization  of  the  volumetric  integral  equations*  In  the  Galerkin 
discretization  scheme  we  expand  the  pressure  in  terms  of  a  set  of  scalar 
locally  supported  trial  basis  functions  0ft, 

N 

p(r)  =  1  o*-3) 

Ct=l 

and  use  the  same  set  as  testing  functions. 

The  details  of  the  derivation  are  described  in  Appendix  B.  Here  we  give 
the  final  result: 

The  discretized  equation  (3.1)  becomes  the  matrix  equation 
N 

^,AQax0  =  ba>  for  a=  1,2,  ,N  ,  (3.4) 

0=1 


3 


where  the  r.h.s. 


ba  =  ^d3r0o(r)p(“c)(r) 

(3.5) 

represents  the  incident  pressure  field,  and  the  elements  of  the  matrix  A  are 

Aa0  =  J  d3r  <f>a{r)  <j>0(r) 

(3.6a) 

k2  J d3rl63r2<l>a{r1)g{r1  r2)  ^1  ^  0fl(r2) 

(3.6b) 

j  d3r1d3r2  [Vri0o(r1)]  g(rl  r2)  •  ^1 

l  [vra^(t 

2)] 

(3.6c) 

=  j  d3r0Q(r)^_  0^(r) 

(3.6d) 

+  k2  /'d3r1dJr20„(rI)p(r1  r2)  (  ^  2'  ,u- 

J  \  «o  p(r2) 

)  Mr2) 

(3.6c) 

|  ^(ra) 

(3.6f) 

Ill 

+ 

%>* 

+ 

%p 

(3.6g) 

In  the  following  we  refer  to  the  terms  AM,  and  ,4°  (given  by  Eqs.  (3.6d) 
-  (3*6f))  as  the  material ,  monopole ,  and  dipole  contributions* 

The  two  alternative  expressions  given  in  Eq.(3,6)  here  are  equivalent. 
The  first  one  is  derived  directly  from  Eq.(3.1),  and  the  second  is  obtained 
by  integrating  by  parts  and  using  the  definition  of  the  Green  function  g  (it 
can  also  be  deduced  from  an  alternative  form  of  the  differential  equation 
given  in  Appendix  A).  However*  the  individual  terms  in  these  equations 
have  quite  different  properties: 

1 .  The  nearly  diagonal  term  (3,6a)  is  independent  of  the  material  prop¬ 
erties,  while  the  term  (3.6d)  (denoted  in  Eq.(3.6g))  docs  depend 
on  p0,  Wc  recall  that  in  our  case  pQ/p  —  10  ~3,  hence  is  small. 

2.  For  a  Large  distance  R  between  the  supports  of  the  basis  functions 
and  both  the  terms  (3.6b)  and  (3.6e)  behave  as  l/R ,  hence  they  can 
be  interpreted  as  due  to  monopole-monopolc  (MM)  couplings.  How¬ 
ever,  for  materials  of  interest,  the  term  (3.6e)  (Aj^)  is  much  smaller 
than  (3,6b)  (typically,  k/kq  ^  10-4). 
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3.  Depending  on  whether  p(r)  is  constant  or  variable  over  the  support 
of  <f>p,  the  term  (3- 6c)  behaves  at  large  distances  as  l/R3  or  l//?2, 
i.e.,  can  be  interpreted  as  a  dipole-dipole  (DD)  or  dipotc-monopolc 
(DM)  coupling.  The  term  (3.6f)  (A^j),  on  the  other  hand,  is  nonzero 
only  for  variable  p(r),  and  then  it  represents  a  dipole-monopole  (DM) 
coupling.  Also  in  this  case  the  term  (3<6f)  (if  nonzero)  is  much  smaller 
than  (3.6c). 

In  other  words,  the  DD  and  DM  couplings  of  the  term  (3.6c)  are  con¬ 
verted  into  MM  and  DM  coupling  of  the  terms  (3,6d)  -  (3.6f).  Generally 
speaking,  the  individual  terms  (3.6a)  -  (3.6c)  tend  to  be  large  compared 
to  the  terms  (3.6d)  -  (3.6f),  which  already  incorporate  substantial  cancella¬ 
tions.  Therefore,  from  the  numerical  point  of  view,  evaluation  of  the  terms 
(3.6a)  “  (3.6c)  requires  a  high  precision,  in  order  to  ensure  sufficient  accu¬ 
racy  after  cancellation.  These  cancellations  arc  already  taken  into  account 
in  Eqs,  (3.  fid)  -  (3.6f),  hence  the  latter  arc  clearly  preferable. 

In  our  present  solver  we  implemented  discretizations  with  piecewise  con¬ 
stant  basis  functions  supported  on  tetrahedra,  and  with  piecewise  linear 
basis  functions  associated  with  vertices  and  supported  on  sets  of  tetrahedra 
sharing  the  given  vertex. 

An  example:  application  to  a  large-contrast  problem.  In  our  appli¬ 
cations  to  acoustic  wave  scattering  on  biological  objects  the  mass/stiffness 
matrix  of  (3.6)  exhibits  a  particular  feature*  the  last  term,  the  DM  compo¬ 
nent  of  (3.6f),  becomes  the  dominant  contribution,  due  to  the  disconti¬ 
nuity  of  the  ratio  p0/p(r)  at  the  interface  of  the  object  and  the  surrounding 
air  (p0/p{r)  =  1  outside  of  the  object  and  p0/p(r)  ~  10“3  inside).  As  we 
discussed  in  the  points  1  and  2  above,  the  remaining  terms,  Aj^  and  Aj^, 
are  small,  because  both  p0/p  and  k/kq  are  small. 

Since  the  surface  term  A^  dominates,  the  system  of  equations  (3.4) 
is  much  more  sensitive  to  the  values  of  the  solution  (the  pressure  p)  on  the 
surface,  than  to  its  values  in  the  interior  of  the  object.  In  other  words, 
the  system  is  ill-conditioned,  and  the  solution  inside  the  object  is  poorly 
defined.  This  fact  is  illustrated  by  the  example  of  a  scattering  on  a  ho¬ 
mogeneous  sphere  with  p/po  =  IQ3  and  k/kq  =  10-4,  and  discretized  with 
N  =  200,000  tetrahedra  (Fig.  1),  and  subject  to  an  incident  pressure  wave 
of  unit  amplitude  and  the  wavelength  A  =  6  m.  In  the  computation  wc  used 
our  fast  solver  code  based  on  the  FFT  compression  of  the  stiffness  matrix, 
described  in  Section  3.3  and  in  Appendix  D, 
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Figure  1:  Distribution  of  the  absolute  value  of  the  pressure!  |p(r)|,  on  a  plane 
parallel  to  the  incident  wave  vector  and  passing  through  the  sphere  center. 
Analytic  solution  is  compared  with  the  conventional  solution  of  the  L-S 
integral  equations  for  p/pu  =  103  and  k/kq  —  10-4,  and  for  a  discretization 
with  N  =  200,000  tetrahedra.  Note  different  pressure  scales  for  the  two 
distributions. 

The  comparison  of  Fig.  1  clearly  shows  that  the  numerical  solution  for 
the  pressure  inside  the  sphere  is  completely  unreliable.  We  note  that  the 
pressure  distributions  in  the  Figure  are  plotted  in  two  different  scales:  up 
to  1.5  for  the  analytic  solution  and  up  to  10  for  the  numerical  solution;  in 
fact,  in  the  latter  case  the  pressure  exceeds  at  some  points  the  value  600. 
By  introducing  small  artificial  “errors”  in  the  components  of  the  r.ihs.  vec¬ 
tor  h  in  Eq.(3.4)  we  also  verified  that  small  fluctuations  in  6  cause  large 
and  unpredictable  changes  in  the  solution  x<  At  the  same  time,  as  wc  have 
have  checked,  the  scattering  cross-section  (which  is,  practically,  sensitive 
only  to  the  pressure  at  the  object  surface)  is  reproduced  with  a  near  perfect 
accuracy.  Both  these  facts  Eire  consistent  with  ill- conditioning  of  the  dis¬ 
cretized  L-S  integral  equations,  resulting  in  a  poorly  determined  solution  in 
the  interior  of  the  object. 

3,2  Development  of  a  novel  two-stage  solution  method  for 
large-contrast  problems 

Reformulation  of  the  integral  equation  system.  We  have  found  that 
the  ill-conditioning  problem  encountered  for  large  density  contrast  can  be 
remedied  by  reformulating,  in  a  rigorous  way,  the  original  integral  equation 
as  a  set  of  a  surface  and  volume  equations,  both  of  them  well  conditioned. 
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The  details  of  our  formulation  are  described  in  Appendix  C;  here  we  only 
sketch  the  solution  procedure: 

We  first  write  the  integral  equation  (3.1)  as 

p(inc)  =  Kp,  (3.7) 

and  define  an  operator  K ^  as  the  p{t)fp  — *  oo  limit  of  the  surface  operator 
By  construction,  this  operator  involves  only  the  values  of  the  pressure  on 
the  object  surface.  Wc  denote  by  p ^  the  solution  of  the  resulting  surface 
equation 

p<inc){r)  =  (/<<°> p<°>)(r)  for  rer«.  (3.8) 

We  can  then  show  (Appendix  C)  that,  if  Eq.(3.8)  holds  on  the  object  surface, 
it  also  holds  for  all  points  r  €  fl,  i.e*  in  the  object  interior. 

Next,  we  represent  the  operator  K  as 

K  =  (3.9) 

with  —  £-1  (K—K^)  and  with  the  parameter  £,  depending  on  the  func¬ 
tion  p(r),  chosen  such  that  the  operator  /f (1)  remains  finite  when  [p(r}|  -*  oo 
(wc  can,  c.g.,  take  £  as  the  average  value  of  |p0/p|  in  the  region  fi). 
Correspondingly,  we  parameterize  the  full  solution  as 

p  =  pI0)+£pF+Pv  .  (31°) 

where  pi0)  is  previously  determined  surface  solution  of  Eq*(3.8),  pi 11  is  a 
correction  to  the  surface  solution  pi°\  and  pv  is  the  volume  part  of  the 
solution,  defined  such  that 

K™pv  =  KsPv  =  0.  (3.11) 

By  substituting  Eqs.  (3.9)  and  (3. 10}  into  Eq*(3.7)  and  making  use  of 
Eqs.  (3.8)  and  (3.11),  wc  find  the  volumetric  equation 

=  —  (-KP<1)rf0>)(r)  for  r  e  H  (3.12) 

for  the  unknown  fields  pa1^  and  pvl  with  the  r.h.s,  expressed  in  terms  of  the 
previously  determined  surface  problem  solution  . 

Finally,  the  entire  solution  procedure  consists  of 

L  solving  the  surface  problem  (3.8), 

2.  using  the  solution  p[0)  to  compute  the  r.h.s.  of  Eq.(3.12), 
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3.  solving  the  volumetric  problem  (3. 12)  for  and  pv;  and 

4.  constructing  the  full  solution  according  to  Eq,(3,10). 

We  stress  that  neither  of  the  equations  (3.8)  and  (3.12)  suffer  from  ill- 
conditioning  in  the  limit  p  — *  oo  (£  — ►  0).  At  the  same  time,  since  no  sniall-£ 
approximations  were  made,  the  procedure  can  be  used  for  any  material  den¬ 
sity.  We  also  emphasise  that  these  equations  are  coupled  in  one  “direction1* 
only,  i.e.,  the  volumetric  equation  depends  on  the  solution  of  the  surface 
equation,  but  not  v.  v>  Hence,  our  system  of  equations  does  not  require  any 
iterative  procedure  involving  alternating  solutions  of  the  surface  and  volume 
problems. 

An  application  to  a  large- contrast  problem.  In  order  to  demonstrate 
the  efficacy  of  our  two-step  solution  procedure,  we  consider  again  the  ex¬ 
ample  discussed  at  the  end  of  the  previous  Section,  3.1  (further  examples  of 
application  of  the  method  are  given  in  Section  3.4).  Now,  by  using  the  two- 
step  solution  scheme,  wc  obtain  a  stable  and  accurate  solution,  as  indicated 
by  Fig.  2. 


Figure  2:  Distribution  of  the  absolute  value  of  the  pressure,  |p(r)|,  on  a  plane 
parallel  to  the  incident  wave  vector  and  passing  through  the  sphere  center, 
for  the  same  material  parameters  and  discretization  as  in  Fig.  1.  As  before, 
the  conventional  and  two-stage  solutions  are  plotted  on  different  scales. 

The  error  distributions  of  the  conventional  and  two-stage  solutions  (com¬ 
pared  to  the  analytical  solution)  arc  plotted  in  the  histograms  of  Fig.  3. 
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Figure  3:  Distribution  of  the  relative  errors  of  the  pressure  on  the  section 
shown  in  Fig.  2,  for  the  conventional  and  two-stage  solutions. 

The  improved  conditioning  of  both  surface  and  volume  problems  is  also 
indicated  by  the  convergence  histories  of  the  iterative  solutions  in  the  con¬ 
ventional  and  two-stage  approach,  shown  in  Fig.  4.  Similar  examples  dis¬ 
cussed  in  Section  3.4  suggests  that  the  conditioning  of  the  matrices  in  the 
two-stage  scheme  is  only  weakly  dependent  on  the  discretization. 
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Figure  4:  Convergence  histories  of  the  iterative  solutions  for  the  conventional 
and  two-stage  schemes,  for  the  high-contrast  sphere  problem  with  AT  = 
200, 000. 
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3*3  Development  of  an  FFT- based  compression  and  fast  so¬ 
lution  methods  for  integral  equations  in  acoustics 

A  crucial  clement  of  our  approach  is  a  fast  iterative  solution  method  uti¬ 
lizing  a  stiffness  matrix  compression  method  based  on  Fast  Fourier  Trans¬ 
forms  (FFTs)  [2].  We  selected  this  technique  rather  than  the  Fast  Multipole 
Method  (FMM)  [1],  since  it  achieves  a  better  performance  in  volumetric  and 
low- frequency  problems,  i.c.,  in  typical  problems  of  acoustic  wave  propaga¬ 
tion  in  biological  media;  it  also  provides  a  smooth  and  uniform  transition  to 
medium  and  high  frequencies. 

We  give  here  a  short  summary  of  the  method;  the  details  are  described 
in  Appendix  D. 

The  FFT*bascd  compression  method  -  the  Adaptive  Integral  Method 
(AIM)  [2]  -  was  originally  developed  in  the  context  of  electromagnetics,  and 
became  one  of  the  most  widely  used  fast  solution  methods.  In  the  present 
effort  wc  modified  it  and  applied  to  general  (variable- density)  acoustics  prob¬ 
lems. 

The  AIM  compression  relics  on  a  decomposition  of  the  stiffness  matrix 
A  of  Eq.(3.4)  into  two  components, 

A  =  ANeai  +  j4Far  (3.13) 

where  the  “ncar-ficld”  term  /lNear  describes  interactions  at  distances  of  sev- 
eral  times  the  spatial  resolution  (or  average  basis  function  support  size),  and 
“far-field”  components  AFat  is  responsible  for  the  remaining  interactions.  By 
construction,  the  ncar-ficld  matrix  ANear  is  sparse,  while  the  (dense)  far- field 
matrix  AFar  can  be  represented  in  a  compressed  representation,  as  a  prod¬ 
uct  of  sparse  matrices  and  a  Tocplitz  matrix  (the  latter  also  becomes  sparse 
upon  taking  its  Fourier  transform,  implemented  as  a  FFT). 

The  compressed  form  of  the  matrix  AF*X  is  obtained  by  representing  the 
original  basis  functions  (say,  i.e+,  the  actual  sources  of  the  field  (here, 
a  pressure  field)  as  “far-ficld-equivalcnt”  distributions  of  point- like  sources 
located  on  nodes  of  a  regular  Cartesian  grid;  the  grid  regularity  gives  rise 
to  the  Toeplitz  property  of  the  resulting  far-ficld  matrix. 

More  specifically,  we  introduce  a  regular  Cartesian  grid  covering  the 
scaitcrcr,  and  characterized  by  some  spacing  h%  comparable  to  the  resolution, 
say  a,  of  the  problem  discretization.  We  then  approximate  the  original  basis 
functions  0Q  by  sets  of  equivalent  point  sources  at  the  nodes  of  an  “expansion 
cube”  Ca)  defined  as  a  set  of  (M  +  1)  x  (M  +  1}  x  (M  +  1)  grid  nodes  u, 
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closest  to  the  center  of  the  support  of  the  basis  function, 

<Mr)  -  0Q(r)  =  X!  “  u>  ■  (3-14) 

u£Ca 

The  set  of  the  (M  +  l)3  expansion  coefficients  V  (strengths  of  the  sources) 
is  determined  so  that  both  functions,  0Q  and  have  the  same  multipole 
moments  up  to  the  order  M.  In  other  words,  the  two  functions  generate  the 
same  (up  to  the  multipole  order  M)  far  fields:  hence  the  concept  of  “far-field- 
equivalent'1  sources, The  criterion  just  mentioned  applies  irrespectively  of  the 
relation  between  the  discretization  resolution  a  and  the  wavelength,  i.c,,  it  is 
valid  for  both  high-frequency  and  subwavelength  problems.  An  alternative 
criterion  based  on  approximate  equality  of  asymptotic  radiated  fields  [4]  may 
be  more  advantageous  at  higher  frequencies,  for  which  a>  A/ 10, 

Fig.  5  shows  schematically  couplings  between  distant  and  near-by  tetra- 
hedra,  giving  rise  to  matrix  elements  and  A^ar,  as  well  as  the  expansion 
cubes  CQ,  and  C 7  associated  with  the  tetrahedra.  It  also  indicates  that 
the  wavelength  A  may  be  much  larger  than  the  discretization  resolution. 

The  far- field  matrix  AFar  is  then  defined  simply  by  replacing  the  basis 
functions  in  A  with  the  basis  functions  (pa .  For  example,  if 

Aa0  -  J  d3rl  d%  <Mrl)  s(rl  -  r2>  <Mr2)  >  (3J5) 


then 

=  j  d3r1d3rs?B(r1).|j{r1~ra)?fl(ra)  s  £  V^fftu-v)  .  (3.16) 

**  u,v 

The  last  expression  shows  that,  indeed,  the  matrix  Apar  is  the  product  of 
two  sparse  matrices  V  and  a  Tocplitz  matrix  g .  The  matrix- vector  multipli¬ 
cations  ANear  x  and  AFar  x  appearing  in  the  iterative  solution  of  the  equation 
(3.4)  can  be  then  implemented  as  sequence  of  sparse- matrix- vector  multipli¬ 
cations  and  of  FFTs,  with  the  overall  storage  O(N)  and  computation  cost 
0(N  log  AT)  (where,  as  before,  N  is  the  total  number  of  unknowns). 
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Figure  5:  A  schematic  representation  of  the  near-  and  far-held  couplings 
between  basis  functions  (tetrahedra)  in  a  model  of  a  human  head.  In  the 
enlarged  view  the  expansion  cubes  are  also  shown.  Typical  discretizations 
are  much  finer  than  shown  here. 

In  our  specific  problem  of  acoustics  with  variable  density  we  use  separate 
expansions  for  the  original  scalar  basis  functions,  for  their  gradients,  and 
for  the  gradients  multiplied  by  material -dependent  factors,  as  appearing  in 
Eqs.  (3.6). 

The  near-field  range  R  and  the  multipole  expansion  order  M  can  be 
adjusted  (together  with  the  grid  spaceing  h)  in  order  to  limit  the  error 
in  the  far-ficld  to  the  desired  level,  and  at  the  same  time  minimize  the 
computational  cost  and  storage.  In  acustics  problems  we  found  that  the 
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optimal  parameters  are 

h  ^  \a  ,  R  cr  3  h  ,  Af  =  2  .  (3.IT) 

These  values  are,  in  fact,  more  favorable  (Lc.,  yield  lower  computational 
costs)  than  in  similar  electromagnetic  problems. 

3*4  Examples  and  numerical  results 

We  present  below  numerical  results  illustrating  the  accuracy,  and  the  capa¬ 
bilities  of  our  code.  The  types  of  problems  are  as  follows: 

•  We  validated  the  code  for  several  problems  involving  homogeneous 
and  layered  spheres  for  the  case  of  constant  density  (heM  the  same 
density  of  the  background  medium  and  the  scatterer)*  We  computed 
numerical  solutions  for  problems  from  few  thousand  to  about  one  mil¬ 
lion  unknowns,  and  found  that  the  numerical  solutions  approach  the 
analytic  one  and  remain  stable  with  the  increasing  discretization  den¬ 
sity. 

•  We  verified  the  accuracy  and  stability  of  our  code  (utilizing,  in  this 
ease,  the  two-stage  solution  scheme  described  in  Section  3.2)  for  a  num¬ 
ber  of  large- contrast ,  variable- density  problems  involving  homo¬ 
geneous  and  layered  spheres  with  material  properties  similar  to  those 
of  biological  tissues. 

•  We  performed  computations  for  a  simple  model  of  a  human  head 
with  the  auditory  channel,  represented  as  a  layered  sphere  with  a  cylin¬ 
drical  channel,  and  composed  of  about  800,000  tetrahedra. 

•  We  carried  out  computations  for  a  realistically  shaped  model  of  a 
human  head  composed  of  approximately  1,000,000  tetrahedra, 

3.4.1  Code  validation  For  constant-density  problems 

In  this  set  of  problems  we  considered  acoustic  media  of  the  same  density 
as  the  surrounding  medium,  i.e.,  p( r)  =  pQ.  We  obtained  numerical  (MoM) 
solutions  by  using  our  fast  FFT  solver. 

We  present  below  results  for  a  three- layer  sphere  of  the  outer  radius  of  L0 
m  and  the  total  number  of  unknowns  of  626,561  and  376,896  corresponding 
to  the  finer  and  to  the  coarser  discretizations  respectively  In  Table  I  we 
list  the  material  parameters  (refraction  index  squared,  n2)  and  the  numbers 
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of  tetrahedra  in  each  layer.  We  stress  that  this  computation  was  carried 
out  only  as  a  test  of  the  solver,  and  is  not  mean  to  represent  any  realistic 
biological  object  in  air  (although  it  could  model  such  an  object  immersed  in 
water). 

The  tetrahedral  mesh  for  the  finer  discretization  is  shown  in  Fig.  6,  The 
sphere  is  excited  by  a  plane  wave  of  wavelength  A  =  2  m. 


Table  1:  Parameters  in  the  three- layer  sphere  problem 


outer  radius 

— IT T~ 
n 

JV(finc) 

N  (coarse) 

0.6  m 

0.20 

153,241 

83,857 

0.8  m 

0.05 

184,029 

109.815 

1.0  m 

0.20 

289,291 

183,324 

total 

626,561 

376,896 

Figure  6:  One  of  the  tetrahedral  nieshes  (with  626,561  tetrahedra)  used  in 
the  computations  for  the  layered  sphere. 

In  Fig.  7  we  plot  the  far- field  (scattering  cross-section)  obtained  using 
the  analytical  solution  and  numerical  solutions  for  the  two  discretizations 
given  in  Table  L  The  curves  are  practically  indistinguishable. 

The  pressure  distribution  inside  the  sphere  is  visualized  in  Fig.  8,  This 
Figure  shows  an  intersection  of  the  tetrahedral  mesh  with  a  plane  passing 
through  the  center  of  the  sphere,  parallel  to  the  incident  wave  vector.  The 
triangular  intersections  of  the  tetrahedra  are  assigned  colors  according  to 
the  pressure  values  in  the  tetrahedron.  The  maximum  difference  between 
the  two  numerical  and  the  exact  analytical  solution  is  less  than  3%. 
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Figure  7:  Comparison  of  the  scattered  field  (cross-section)  for  the  exact 
solution  for  a  layered  sphere  and  our  numerical  MoM  solutions  for  N  = 
626,561  and  N  =  376,896  unknowns. 
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Figure  8:  Comparison  of  pressure  distributions  (the  imaginary  part  of  p(r ) ) 
in  a  layered  sphere,  computed  using  the  exact  and  numerical  MoM  solutions. 
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3.4.2  Code  validation  for  variable-density,  large-contrast  prob¬ 
lems 

Any  problem  involving  scattering  on  acoustic  objects  in  air  must  necessarily 
involve  large  ratios  of  densities.  In  particular,  the  ratio  of  density  of  a 
biological  tissue  (or  water)  to  that  of  air  is  above  1000.  Accordingly,  in  our 
validation  tests  wc  considered  such  a  density  and  the  refraction  index  in 
the  range  0  <  n  <  1,  also  typical  of  biological  media  (at  least  for  pressure 
waves). 

Homogeneous  spheres.  In  order  to  verify  the  accuracy  of  our  acoustic 
integral-equations  solver,  wc  carried  out  a  number  of  computations  for  either 
homogeneous  or  layered  spheres,  and  compared  the  results  with  the  exact 
analytical  solutions,  expressed  as  series  involving  Bessel  functions  and  Leg¬ 
endre  polynomials.  The  use  of  the  developed  two-stage  solution  scheme 
was  crucial  in  obtaining  the  results,  otherwise  the  problems  would  become 
ill-conditioned  and  the  solutions  unstable. 

Figs,  9  to  10  show  results  for  a  sphere  of  radius  a  ~  10cm,  immersed  in 
air,  and  filled  with  an  acoustic  material  of  the  relative  density  p/pQ  =  1G3 
and  relative  compressibility  k/k0  —  IQ-4,  corresponding  to  n 2  =  0.1  or  n  ~ 
0.3 16.  These  parameters  are  similar  to  those  of  bone  or  water  (sec  Table  3 
of  Section  3.5.2).  In  the  numerical  solution  we  use  the  sphere  discretization 
with  about  N  =  970, 000  tetrahedra.  The  sphere  is  excited  by  a  plane  wave 
of  wavelength  A  —  60cm  (he,,  of  frequency  /  ~  573 Hz). 

The  current  version  of  the  solver,  partly  parallelized  for  shared-memory 
(MF)  platforms,  was  run  on  an  SGI  Origin  3900  system.  The  total  compu¬ 
tation  time  on  8  processors  was  about  130  minutes.  The  maximum  required 
memory  was  about  6.4  GB. 

In  Fig.  9  we  plot  distributions  of  the  absolute  value  of  the  pressure  on  a 
plane  parallel  to  the  incident  wave  vector.  Fig.  10  shows  the  histogram  of 
the  relative  error  in  the  pressure  of  the  numerical  solution  versus  the  ana¬ 
lytical  one.  The  histogram  is  based  on  31406  events  representing  pressures 
at  polygons  obtained  by  intersecting  the  tetrahedral  mesh  with  the  plane  of 
Fig.  9.  It  is  seen  that  the  error  in  the  two-stage  solution  is  of  the  order  of 
1%. 

In  Fig,  10  wc  compare  the  stage-1  and  stage-2  convergence  of  the  two- 
stage  iterative  solution  scheme  with  that  of  the  conventional  scheme.  The 
slow  convergence  of  the  conventional  procedure  (and  the  resulting  large  fluc¬ 
tuations  in  the  pressure  distributions  (not  shown  here))  are  typical  symp¬ 
toms  of  poor  conditioning.  These  difficulties  do  not  appear  in  the  two-stage 
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solution  approach.  Incidentally,  we  note  here  that,  while  the  convergence  of 
conventional  solution  in  the  present  problem  is  significant  slower  than  for 
the  smaller  problem  with  200.000  unknowns  (Fig.  4),  the  numbers  of  iter¬ 
ations  in  the  two-stage  scheme  remain  almost  unchanged,  which  indicates 
that  the  conditioning  of  the  surface  and  volume  operators  in  Eqs.  (3.8)  and 
(3.12)  is  largely  independent  of  the  discretization, 

analytical  solution  numerical  solution 

•  i  •  > 


20  cm  3  X  /  3 

Figure  9:  Distribution  of  the  absolute  value  of  the  pressure,  |p(r)|,  on  a 
plane  parallel  to  the  incident  wave  vector,  passing  through  the  homogeneous 
sphere  center,  computed  analytically  and  by  means  of  our  acoustic  solver, 
for  a  discretization  with  N  —  970, 000  tetrahedra. 


Figure  10:  Histogram  of  the  relative  error  in  the  pressure  for  the  section 
shown  in  Fig,  9, 
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Figure  11:  Convergence  of  the  two  iterative  solutions  in  the  two-stage 
scheme,  compared  to  the  convergence  of  the  the  conventional  solution. 


Multi-layered  spheres*  In  Figs.  12  and  13  we  show  results  for  the  pres¬ 
sure  distribution  in  a  layered  sphere,  with  the  outer  shell  chosen  to  represent 
the  bone,  and  the  interior  of  the  sphere  described  by  mechanical  parame¬ 
ters  of  water  (Tables  2  and  3).  According  to  Table  3,  the  parameters  of 
the  inner  region  approximate  the  behavior  of  pressure  waves  in  the  brain 
tissue.  As  shown  in  Fig.  12,  the  outer  radius  of  the  sphere  is  10  cm,  and  the 
inner  radius  9  cm.  We  constructed  a  tetrahedral  mesh  with  the  tetrahedron 
sizes  about  4  mm;  hence  the  thickness  of  the  shell  is  about  three  times  the 
tetrahedron  size.  The  total  number  of  tetrahedra  is  about  N  ^  755,000. 


Figure  12:  Discretization  of  the  layered  sphere  with  N  =  755, 000  tetrahedra. 


Fig,  13  shows  distribution  of  the  absolute  value  of  the  pressure,  |p(r)  for 
the  incident  wave  with  A  —  60  cm  (frequency  /  ^  573  Hz),  As  in  the  case 
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of  the  homogeneous  sphere,  the  errors  in  the  computed  pressure  arc  of  the 
order  of  1  %* 


Figure  13:  Distribution  of  the  absolute  value  of  the  pressure,  |p(r)|,  on 
a  plane  parallel  passing  through  the  layered  sphere  center,  computed  an¬ 
alytically  and  by  means  of  our  acoustic  solver,  for  a  discretization  with 
N  =  755,000  tetrahedra. 

3*4,3  Other  variable  density  problems* 

We  present  here  some  other  representative  examples  illustrating  capabilities 
of  our  solver  in  problems  characterized  by  large  density  contrast. 

Results  for  a  layered  sphere  with  a  cylindrical  channel*  We  also 
carried  out  computations  for  a  layered  sphere  with  a  cylindrical  channel 
-  a  simple  model  of  an  car  canal.  The  sphere  size,  material  parameters 
of  the  layers,  and  the  frequency  are  the  same  as  in  the  previous  example. 
The  object  is  discretized  with  about  N  =  763,000  tetrahedra,  as  shown  in 
Fig,  14,  The  obtained  distribution  of  pressure  is  plotted  in  Fig.  15,  The 
computational  cost  and  the  convergence  of  the  solution  was  practically  the 
same  as  in  the  previous  examples* 


Figure  14:  Discretization  of  the  layered  sphere  with  a  channel,  with  N  = 
763*000  tetrahedra. 
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1.5 


1.4 


Figure  15:  Distribution  of  the  absolute  value  of  the  pressure,  |p(r)|,  on  a 
plane  parallel  passing  through  the  object  center,  computed  with  the  dis¬ 
cretization  of  Fig.  14. 

Computations  for  a  head  model*  As  an  example  of  a  more  realistic 
geometry,  we  considered  a  human  head  model,  represented  originally  as  a 
surface  mesh  with  about  43,000  triangular  facets  (Fig.  16),  and  constructed 
a  tetrahedral  mesh  with  about  N  =  1,090,000  tetrahedra* 

We  took  the  same  material  parameters  and  the  same  incident  wave  fre¬ 
quency  as  in  the  previous  sphere  problem  (he,,  p/pQ  =  103,  k./kq  ~  IQ-4, 
n  0.316,  A  =  60cm,  /  ~  573 Hz).  At  this  frequency,  the  average  edge 
length  is  about  A/200, 

Fig.  17  shows  the  distribution  of  the  absolute  value  of  the  pressure,  |p(r)|, 
in  the  axial  plane,  for  the  indicated  direction  of  the  incident  wave.  As  before, 
the  solution  was  obtained  by  means  of  the  two-stage  scheme.  Convergence 
of  the  solutions  was  similar  to  that  in  the  sphere  problem. 
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Figure  16:  The  surface  model  of  the  human  head  and  its  representation  as 
a  tetrahedral  mesh  (with  N  =  1,090,000  tetrahedra),  shown  in  the  sagittal 
plane  cut. 


20  cm  =  X  /  3 


Figure  17:  Distribution  of  the  absolute  value  of  the  pressure,  ]p(r)|,  in  the 
axial  plane,  for  the  head  model  of  Fig.  16.  The  computation  was  done  for 
N  —  1,090,000  unknowns  with  the  two-stage  solution  scheme.  A  section  of 
the  tetrahedral  mesh  is  shown  to  indicate  the  spatial  resolution. 
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3.5  A  need  for  elasticity  modeling 

Wc  conclude  the  review  of  our  Phase  I  results  with  a  discussion  on  the  range 
of  applicability  of  acoustic  modeling  of  biological  tissues,  and  the  possible 
need  of  clastic  modeling. 

3.5.1  Properties  of  biological  tissues 

There  exists  a  large  body  of  experimental  data  related  to  biomechanics  of 
the  human  head,  collected  in  various  contexts:  (ultra)sound  propagation, 
impact  mechanics,  neurosurgery,  etc.  (e.g,,  [5,  6]).  We  summarize  here  the 
most  important  mechanical  (acoustic  and  elastic)  properties  of  the  relevant 
tissues. 

Skull  bones.  It  is  generally  accepted  that  skull  bones  [7]  should  be  mod¬ 
eled  as  clastic  media.  It  is  also  known  that  the  bone  properties  may  sub¬ 
stantially  vary  between  individuals,  and  depend  on  the  location  within  the 
given  skull.  Similarly,  soft  tissues  arc  usually  described  as  elastic  materials, 
but  may  be  approximated  as  acoustic  media. 

Brain  tissue.  A  mode  difficult  (and  somewhat  controversial)  issue  is  that 
of  the  constitutive  model  of  the  brain  tissues  (sec  the  review-  [6]).  Here, 
within  the  scope  of  linear  materials,  solid  or  fluid,  elastic,  viscoelastic,  or 
poroelastic  models  are  being  used.  These  models  differ  mainly  in  their  long¬ 
time  behavior  (presence  or  absence  of  shape  memory);  this  differences,  how¬ 
ever,  arc  not  critical  in  our  problem,  in  which  short-time  (nearly  impulsive) 
excitations  arc  involved,  and  wc  do  not  have  to  describe  large-time  scale 
behavior  of  the  system.  On  the  other  hand,  correct  modeling  of  viscous 
properties  of  the  brain  tissue  is  important  in  order  to  describe  energy  dissi¬ 
pation  and  thus  wave  attenuation  in  the  tissue, 

A  related  question  is  that  of  the  brain  tissue  compressibility.  Since  the 
brain  is  about  83%  water  by  weight,  it  should  be  expected  to  be  nearly 
incompressible.  According  to  the  standard  set  of  data  (Refs.  [8,  9,  10,  11]), 
the  bulk  modulus  of  the  tissue  is  K  ^  2*  10y  Pa,  about  IQ5  times  larger  than 
its  shear  modulus  G,  hence  the  Poisson  ratio  u  is  very  close  to  ^  (see  the 
summary  of  mechanical  material  parameters  below-,  in  Section  3.5.2). 

An  additional  problem  is  a  possible  relevance  of  non-linear  effects:  be¬ 
cause  to  the  low  shear  wave  speed  in  the  brain  tissue,  material  velocities 
may  become  comparable  to  phase  (wave- front)  velocity,  resulting  in  steep¬ 
ening  of  wave- fronts,  wave  overturning,  and  similar  phenomena.  However, 
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while  these  mechanisms  may  be  relevant  in  the  case  of  large  impact  forces 
and  Large  displacements  (potentially  causing  head  injury) ,  they  do  not  seem 
important  in  the  context  of  sound  conduction  at  amplitude  levels  that  do 
not  cause  permanent  damage  to  the  tissues. 


Sub-arachnoidal  space.  An  important  element  in  the  human  head  me¬ 
chanics  is  the  sub-arachnoidal  space  between  the  skull  and  the  brain,  filled 
with  the  spongy  trabeculae  and  cerebrospinal  fluid,  and  cushioning  the  brain 
from  shocks  and  vibrations.  This  tissue  is  usually  modeled  as  a  viscous  solid 
of  low  shear  modulus. 

The  properties  of  the  tissues  listed  above  indicate  that  the  entire  human 
head  can  be  described  by  means  of  visco-elasticity  equations.  In  the  follow¬ 
ing  we  argue  that  visco-elastic  modeling  may  be  also  necessary,  i.c.,  a  purely 
acoustic  model  may  not  reproduce  important  properties  of  the  system. 


3.5,2  Elastic  vs,  acoustic  modeling  of  tissues 

We  give  here  a  brief  account  of  the  mechanical  parameters  of  the  pertinent 
materials  -  the  tissues  constituting  the  human  head.  Instead  of  the  Lame 
coefficients  AL  and  f.iL,  we  use  here  the  bulk  and  shear  moduli  K  and  G, 

K  =  K+2H  =  -J^iE-y  0  =  „l  =  _A_,  (3.18) 

in  terms  of  which  the  Young  modulus  and  the  Poisson  ratio  arc  given  by 


E  = 


(3 K  -  4G)  G  _  (3cp  -  4c|)  c§ 


K  -  G 


r2  _  r'2 
Lp  Lg 


and 


v  = 


K  -  2 G  _  4  - 


2 (K-G)  2(4-4)  ‘ 

The  corresponding  speeds  of  sound  arc  then 


Cp  — 


(l-v)E 


(1  + 1/)  (1  -  2v)  p 


tuid 


2(l  +  i t)p  ' 


(3.19a) 


(3.19h) 


(3.20a) 


(3.20b) 
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where  the  subscripts  P  and  S  refer  to  the  pressure  and  shear  waves. 

The  Tables  2  and  3  below  give  the  approximate  mechanical  parameters 
for  some  relevant  materials.  To  the  first  approximation,  soft  tissues  of  the 
head  (skin,  muscles,  etc.}  have  properties  similar  to  those  of  water;  however, 
properties  of  the  brain  tissue  are  quite  different. 

In  the  case  of  the  brain  we  show  both  the  real  parts  of  the  material 
parameters  (in  the  first  line),  and  imaginary  parts  (in  the  second  line);  for 
other  tissues  the  imaginary  parts  of  the  parameters  are  much  smaller.  It 
is  known  that  the  large  imaginary  parts  [9]  arc  important  in  describing 
viscoelastic  properties  of  that  tissue. 

In  Tabic  3  we  list  indices  of  refraction,  n?  -  %/cp  and  ns  =  e0/cs, 
relative  to  the  environment  medium  (air),  for  which  Cq  —  344  m/s.  The 
effective  acoustical  compressibilities  and  /cs,  arc  defined  such  that 

pKp  =  Cp  ,  pKs  —  c|  .  (3.21) 

Their  values,  relative  to  the  compressibility  of  air,  as  well  as  relative  medium 
densities,  are  also  listed;  again,  for  the  brain  tissue,  we  show  imaginary 
parts  of  the  parameters  below  the  real  parts.2  The  parameters  of  Table  3 
are  useful  in  comparing  elasticity  problems  with  approximately  equivalent 
acoustics  problems. 

Table  2:  Mechanical  properties  of  some  materials 


material 

P  !g/m3I 

B  [Pa] 

V 

K  [Pa] 

G  [Pa] 

air 

1.2 

1.42  e+05 

0 

water 

1000 

0.5 

2.25  e+10 

0 

bone 

2132.6 

1.379  c+09 

0.25 

1.65  e+09 

5.516  e+08 

brain 

Im: 

1002.5 

1.752  e+05 
—3.630  e+05 

0.49999 
0.39  e-04 

2.103  c+09 
-1.613  e+05 

5.840  c+04 
-1.210  e+05 

rubber 

1100,0 

8.0  e+06 

0.333 

1.2  c+06 

3.0  c+05 

2Since  the  brain  tissue  is  dispersive,  these  quantities  also  depend  on  the  frequency.  The 
listed  values  correspond  to  1  kHz. 
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Tabic  3:  Speeds  of  sound,  relative  refraction  indices,  and  equivalent  com¬ 
pressibilities  of  some  materials 


material 

cp  [m/s] 

cs  [m/s] 

nP 

ns 

p/p o 

KY>fKQ 

k$/kq 

air 

344 

1 

water 

1500 

0.22993 

833.33 

6.311  e-05 

bone 

880.88 

508.58 

0.3905 

0.6764 

1777.2 

8.5813  c-05 

2.574  c-04 

brain 

!m: 

1448.4 

-0.0556 

9.804 

-6.155 

0.2375 
9.1  c-06 

25.167 

15.799 

835.42 

835.42 

6.752  e-05 
5. 179  c- 09 

0.4594 

0.9518 

rubber 

33.03 

16.514 

10.41 

20.83 

916.67 

0.1183 

0.4733 

In  addition  to  the  material  parameters  of  the  tissues  we  included  typical 
parameters  of  rubber,  in  order  to  contrast  them  with  the  those  of  the  brain. 
Although  both  materials  are  highly  elastic,  the  difference  is  that  the  brain 
tissue  is  nearly  incompressible  (hence  a  large  bulk  modulus  K  (Table  2), 
much  larger  than  the  shear  modulus  G ),  while  rubber  is  much  more  com¬ 
pressible  and  its  bulk  modulus  is  only  about  5  times  larger  than  the  shear 
modulus.  As  a  consequence,  the  difference  between  the  pressure  and  shear 
wave  speeds  in  brain  is  much  greater  than  in  rubber  (Table  3), 

Another  important  aspect  of  the  mechanical  properties  of  the  brain  tis¬ 
sue  is  its  viscosity,  manifesting  itself  mostly  in  the  large  imaginary  part  of 
the  shear  modulus  (Table  2)  and,  consequently,  in  the  shear  wave  speed 
(Tabic  3). 

3.5,3  Range  of  applicability  of  acoustic  modeling  of  tissues 

We  have  considered  the  question  of  what  aspects  of  the  problem  of  sound 
conduction  in  the  human  head  could  be  described  in  terms  of  acoustics t 
without  invoking  clastomechanics  at  all  Our  assessment  is  based  on  the 
observation  that  the  presence  of  slow  shear  waves  in  the  brain  (and  thus 
description  of  the  brain  in  terms  of  (visco)  elasticity)  appears  to  be  necessary 
in  order  to  reproduce  the  experimentally  measured  resonance  and  antircso- 
nance  structure  of  the  pressure  distribution  as  a  function  of  frequency. 

One  of  the  reasons  which  suggest  importance  of  shear  waves  in  tissues  is 
that,  if  only  fast  pressure  waves  (cp  >  1500  m/s)  arc  taken  into  account,  the 
lowest  resonant  structure  in  the  human  skull  filled  with  the  brain  or  other 
sort  tissue  appears  to  arise  only  at  about  6  kHz,  while  the  lowest  resonances 
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in  the  living  skulls  arc  observed  in  the  vicinity  of  1  kHz  [12].  This  fact 
suggests  that  the  measured  resonance  structure  is  associated  with  waves  of 
much  lower  speeds  and  thus  much  smaller  wavelengths. 

To  be  more  precise,  resonances  below  1  kHz  were  observed  also  for 
empty  dry  skulls  (or  skulls  artificially  covered  with  a  relatively  thin  layer 
of  a  damping  substance  [13]}.  However,  if  the  shell  is  filled  with  a  low- 
compressibility  acoustic  material  (such  as  water),  the  system  becomes  ef¬ 
fectively  so  “stiff”3  that  the  lowest  resonance  appears  only  at  about.  7  kHz. 
This  result  suggests  that,  in  order  to  generate  the  observed  resonance  struc¬ 
ture,  the  material  filling  the  shell  (the  skull)  must  support  slower  waves.  In 
reality,  such  waves  could  be  the  shear  waves  supported  by  the  brain  tissue, 
of  speed  less  than  10 m/s  (Table  3). 

To  summarize,  our  assessment  is  that 

■  Purely  acoustic  modeling  may  be  reasonable  for  (visco)clastic  tissues 
in  which  the  pressure  wave  and  shear  wave  speeds  arc  comparable,  i.eM 
the  the  Poisson  ratio  v  is  not  too  close  to 

■  However,  tissues  in  which  the  Poisson  ratio  is  close  to  and  in  which 
the  shear  waves  are  much  slower  than  the  pressure  waves,  require  full 
clastic  modeling;  it  does  not  seem  that  modeling  of  shear  waves  as 
acoustic  pressure  waves  can  give  correct  results. 

Consequently,  we  believe  that  the  appropriate  approach  should  be  based 
on  elasticity,  or,  for  the  brain  tissues,  viscoelasticity,  which  should  include 
energy  dissipation  effects* 


3i.e.,  the  speed  of  sound  waves  becomes  high. 
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Appendices 


A  Integral  equations  in  acoustics 


The  starting  point  in  deriving  the  integral  Lippmaim-Schwinger  equations 
in  acoustics  is  the  frequency-domain  differential  equation  [3] 

w2«(r)p{r)  +  V  ■  Vp(r)j  =0  ,  (A.l) 

in  which  u  is  the  frequency,  p  the  equilibrium  density  of  the  medium,  and 
K  ist  compressibility,  and  p  the  acoustic  (excess)  pressure. 

The  L  ipp  mann-  S chwinger  (L-S)  equation  for  the  pressure  can  be  ob¬ 
tained  in  the  standard  way  from  the  differential  equation  (AT)  rewritten  in 
the  form 


(, fc 2  +  V2)  p  -  k2  ^1  -  p  _  V  ■  Vp  =0  (A. 2) 


with 


k  =  —  =  VPo«o 
Co 


(A.3) 


being  the  wave  number  of  the  incident  wave  in  the  background  medium. 
Eq.(A.2)  leads  [3]  to  the  L-S  equation  (Eq.(3.1)  in  the  text) 


n(mc)(r)  _  +  J  d3f»  |/;2s^r  _  r'j  j^i  _  Mil  j  p(r') 

-[Vr,S(r-r')]-(l-^y)  Vr,p(r')j, 


(A.4) 


where  p("ic)  is  the  incident  wave,  satisfying  the  Helmholtz  equation  in  the 
background  medium, 


(k2  +  V2)  p<inc!(r}  =  0 
and  g  is  the  corresponding  Green  function, 

eifc|r| 


S(r)  = 


4?r|r| 


(A.5) 


(A.6) 
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Derivation  of  Eq.(A*4)  involves  the  definition  (A. 5}  of  the  incident  wave, 
the  standard  representation  theorem  for  the  field  p,  and  the  Sommerfcld 
radiation  condition  imposed  on  the  scattered  wave  p^sc*  =  p  —  p(incK 

In  the  main  text  wc  gave,  in  Eq.(3.6),  two  alternative  expressions  for 
the  matrix  elements  of  the  stiffness  matrix*  The  first  of  these  was  derived 
directly  from  Eq,(A*4).  The  second  can  derived  by  integrating  by  parts,  or 
can  foe  obtained  from  an  alternative  L-S  equation 

pp.0,(r)  _  J?Lp(r)  +  y  d V  jfc2s(r  -  r')  P(r') 

^  v  (A.7) 

-[Vr.9(r-r')].(vr,^jKr')|. 

which  follows  from  the  differential  equation  (A, 2)  rewritten  in  the  form 


(t*  +  v2)  (£a„)  +  *2  (i-a)  p  +  v 


(A.8) 


B  Discretization  of  integral  equations  in  acoustics 

B.l  Discretization  with  piecewise  constant  basis  functions 

One  of  the  simplest  choices  of  the  basis  functions  in  Eq,(3,3)  are  constant 
functions  4>a  (of  value,  sayT  X)  supported  on  tetrahedra  ta.  Wc  also  assume 
that  the  material  parameters  p(r)  and  rc(r)  arc  constant  within  individual 
tetrahedra.  Under  these  assumptions  the  matrix  component  AR  is  diagonal, 
and  its  elements  (3.6d)  are  proportional  to  the  volumes  of  the  tetrahedra. 
The  matrix  elements  (3*6c)  arc  given  by  integrals  over  tetrahedra,  and  the 
matrix  elements  (3.6f)  are  expressed  as  sums  of  integrals  over  faces  of  the 
tetrahedra,  since  gradients  of  the  basis  functions  arc  proportional  to  delta- 
functions  on  the  tetrahedron  boundaries.  Explicit  expressions  arc 


a 

^.1  a 

c 

* 

II 

(B.la) 

/  A  j 

^%s(r,-r2)  , 

(B.  lb) 

and 
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a°/3=  J  dS  j  d2r2j(ri  -  r2)nQ(r1)-nja(r2)A(r2) 

«o  Mg 

=  53  5Z  ft/u ‘”/3  A/,  /tl2ri  /d2r2  5(ri  ~r2)  ■  (B-lc) 
/Qesea/fl6*a  “  J/n  /fl 


In  these  expressions  is  the  volume  of  the  tetrahedron  ia;  pa  and  /ca  denote 
values  of  these  parameters  on  the  tetrahedron  ta\  hQ(r)  is  the  outer  unit 
normal  to  the  the  tetrahedron  tQ\  fQ  is  a  face  of  the  tetrahedron,  and  n 

fa 

is  the  outward-oriented  normal  to  that  face;  finally,  A  is  the  discontinuity 

/  a 

of  the  function  p^/pir)  across  the  face 


lim  | 
£■ — *0  + 


Pa 


Pa 


A  r  +  enfJ  p(  r-fn;J, 


for  r  €fa 


(B.2) 


The  sums  over  faces  in  Eq.(B.lc)  arise  because  gradients  of  the  basis  func¬ 
tions  arc  proportional  to  delta-functions  on  the  tetrahedron  boundaries.  Ac¬ 
tually,  we  were  also  able  to  convert  volume  integrals  to  non-singular  surface 
integrals. 

Clearly,  the  matrix  element  is  nonzero  only  if  the  tetrahedron  t(j  is 
adjacent  to  a  surface  S  across  which  p(r)  is  discontinuous.  In  particular,  in 
the  case  of  a  single  volume  fl  characterized  by  a  constant  density  p,  immersed 
in  the  background  medium,  Eq.(B.2)  yields 


for  boundary  face  /  , 
for  interior  face  /  . 


(B.3) 


We  have  implemented  in  our  code  discretization  based  on  the  constant 
basis  functions  described  here  and  found  that,  in  spite  of  their  simplicity, 
accurate  solutions  of  large  problems  can  be  obtained. 


B.2  Discretization  with  piecewise  linear  basis  functions. 

As  another  choice  of  basis  functions  in  acoustic  L-S  equations  we  considered 
scalar  piecewise  linear  functions  associated  with  a  vertex  of  the  tetrahedral 
mesh  and  supported  on  the  set  of  tetrahedra  sharing  that  vertex  (Fig,  IS). 
It  is  uniquely  defined  by  linearity  and  by  the  requirement  that  0tf(rv)  =  1 
and  <£v(r)  =  0  for  r  on  any  external  face  of  any  of  the  tetrahedra.  It  is  also 
continuous  on  the  region  defined  by  Tv. 
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Figure  18:  The  set  of  tetrahedra  supporting  the  basis  function  <j>v  associated 
with  the  vertex  v.  A  tetrahedron  t  (one  of  the  set  7\n  consisting,  in  this  ease, 
of  10  tetrahedra)  is  shown  highlighted,  together  with  the  “height  vector’5  ht , 
normal  to  the  face  (vi,«2i03)< 

In  the  following  we  identify  the  index  a  in  <j>a  with  the  appropriate 
vertex  number-  In  this  formulation  we  also  retain  the  previous  assumption 
of  material  parameters  constant  on  tetrahedra. 

The  matrix  elements  A^,  A^,  and  A^  can  be  now  evaluated  in  a 
similar  way  as  for  constant  basis  functions.  One  of  the  differences  is  that  the 
material  term  of  the  matrix,  AR,  is  no  longer  diagonal,  but  remains  nearly 
diagonal.  The  gradient  of  appearing  in  A^  (Eq.(3.6f))  is  nonsingular, 
since  the  basis  functions  themselves  arc  continuous.  It  can  also  be  shown 
that  tetrahedra  adjacent  to  interfaces  between  different  materials  do  not 
require  any  special  treatment. 

The  advantage  of  the  linear  basis  functions  is  that  they  provide  a  smoother 
solution  for  pressure  and  significantly  reduce  the  number  of  unknowns 
in  the  problem,  compared  to  tetrahedron-based  functions.  The  reason  is 
that,  for  typical  tetrahedral  meshes  the  number  of  vertices  is  much  smaller 
than  the  number  of  tetrahedra,  typically 


C  A  two-stage  solution  method 

In  utilizing  our  discretized  L-S  acoustics  integral  equations  we  encountered 
a  serious  difficulty  associated  with  poor  conditioning  of  the  resulting  linear 
system:  As  we  discussed  in  Section  3.1  and  in  Appendix  B,  a  peculiar  feature 
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of  the  integral  equation  (A* 7)  (and  thus  the  resulting  equation  involving  the 
matrix  dements  (3.6)}  is  that,  for  |p(r)/p0|  1  and  moderate  values  of 

the  refraction  index  n(r),  the  last  (surface)  term,  resulting  in  the  dipole 
term  AD  in  the  matrix,  dominates  over  the  remaining  terms.  It  follows 
that  the  matrix  equation  becomes  ill-conditioned:  the  surface  component  of 
the  solution  is  well  defined,  but  the  volume  part  of  the  solution  (pressure 
in  the  interior  of  the  object)  is  poorly  determined.  This  is  what  we  found 
numerically  in  applying  the  initial  version  of  our  solver  to  problems  with 
large  material  density  contrasts. 

We  were,  however,  able  to  remedy  the  difficulty  and  turn  the  situation  to 
our  advantage  by  separating  the  solution  into  its  surface  and  volume  parts, 
renormalizing  them,  and  solving  the  problem  in  two  stages,  both  without 
introducing  any  approximations  (such  as  an  expansion  in  pQ/p). 


CM  Reformulation  of  integral  equations 

In  order  to  separate  the  surface  and  volume  contributions  to  the  solutions, 
we  represent  the  integral  equation  (A. 7)  as 

p(mc)  _  p  =  j Ksp  +  A'v p  ,  (CM } 


where  the  surface  and  volume  operators  arc  defined  by 

(Ksp)(r)  =  ~  f  d V  [n(r')  -  Vr,g( r  -  r')]  A(r')?>(r') 
Jan 


Ian 


(C.2a) 


and 


//’■'«<*  -  o  -  $j]  pm 

<c-2b) 


The  surface  operator  Ks  of  Eq.(C,2b)  is  defined  as  due  to  the  dominant 
discontinuity  of  the  density  on  the  object  boundary  9fi;  the  volume  operator 
Kv  includes  contributions  due  to  variations  of  density  in  the  interior  of  the 
object. 

In  the  following  we  will  obtain  a  set  of  integral  equations  equivalent 
to  Eq.(C.l),  but  characterized  by  better  conditioning  in  the  limit  of  high 
object  density,  p/pa  S>  1.  In  the  derivation  we  will  utilize  the  standard 
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representation  theorem  for  exterior  scattering  problems  in  acoustics, 


for  r  6  R3  \  0  , 
r)  for  r  6  , 

for  r  €  ft  . 

(C.3) 

Eq.  (C.3)  is  valid  (under  appropriate  regularity  conditions  on  the  surface 
30,  see,  e.g.,  Ref.  [14])  for  any  function  p  (here  interpreted  as  the  total 
pressure  field)  satisfying  the  Helmholtz  equation  in  R3  \0,  and  approaching 
p*m^  at  infinity,  in  the  sense  that  the  scattered  field  =  p-p(inc)  satisfies 
the  Sommerfeld  boundary  condition, 

f-Vp^(r)  -  ikj^(r)  =  o(jrj“l)  for  |r|  -*  oo  .  (C.4) 

Eq,(C-3)  follows,  in  fact,  from  applying  the  second  Greco  identity  to  the 
functions  g  and  p  —  p^incK 

Our  first  step  in  reformulating  the  integral  equation  (C.l)  is  to  define 
an  operator  as  the  p(r)/p  — *  oo  (i.e.,  A(r)  — *  1)  limit  of  the  surface 
operator  Ksi 

(KW  p)(r)  =  -  [  dV  3g„(r~r'}  p(r')  for  r  6  n  ;  (C.5) 

Jan  ®n(T ) 

the  operator  K ^  maps  functions  p  supported  on  dSl  onto  functions  defined 
on  fi,  In  terms  of  this  operator,  we  then  define,  on  the  boundary  Ofi,  a 
function  psU)(r)  as  the  solution  of  the  surface  integral  equation 

pfinc)(r)  -  {K^p^)(r)  for  r  G  3fi_  ,  (C.6) 

where  dSl_  C  fi  is  a  surface  infinitesimally  close  to  the  boundary  3fi  and 
located  inside  the  region  fi. 

At  the  same  time,  if  p  is  the  solution  of  the  exterior  Neumann  problem 
with  the  hard'Surface  boundary  condition 

,  r)  =0  for  r  €  an  ,  (C.7) 

3n(r) 

and  with  the  incident  w^avc  p^nc\  the  representation  theorem  (C.3)  yields 
p(incV)  +  /  d2r'  ^g(r— ^  p(r')  =  0  for  r  e  fi  ,  (C.8) 

Jdtt  dn(v ) 


p(mc)(r)+  f  ctar'  p(r’)-n(r-r‘) 

Jon  L  3n(r) 


dp{r') 


dn{  r1) 


'p(r) 

0 
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Eq.(C-S)  is  thus  equivalent  to  Eq+{C6),  provided  wo  identify  with  p  on 
It  also  follows  that,  if  Eq^(C.G)  holds  on  the  surface  dfl_ ,  it  also  holds 
everywhere  in  the  region  ft,  i.e.,  Eq.(C.6)  implies 

pOnc)(r)  =  (tfW)  p£°>)  (r)  for  r  €  ft  .  (C.9) 

Our  second  step  in  the  solution  procedure  is  to  represent  the  operator 
A'  as 

K  =  K{0)+£K{1)  (C.10) 

with 

!<w  =  £-' (K  -  1<10))  ,  (C.ll) 

and  with  the  parameter  f,  depending  on  the  function  p(r),  chosen  such  that 
the  operator  remains  finite  when  |p(r)|  — *  oo.  More  explicitly 


p)(r)  =  f  dV 

Jm 


dg{  r-r1) 
tJn(r') 


Po 

£p(r') 


p(r') 


and  wc  can  choose  c.g.,  as  the  average  value  of  pa/p{ r)  over  the  region  fi, 

•  (C.12) 

Having  specified  the  operators  and  K^\  and  the  solution  p^\  wc 
represent  the  full  solution  as 

P  =  pi0)  +€p*1j  +Pv  .  (C.13) 

where  pi’1  is  the  correction  to  the  surface  solution  p£0)  and  pv  is  the  volume 
part  of  the  solution,  defined  such  that 

Pv  =  KgPy  =  0  .  (C.14) 

By  substituting  Eqs.  (C.10)  and  (C.13)  into  Eq.(C.l)  and  making  use  of 
Eqs.  (C.G),  (C.9),  and  (C.14),  we  find  the  volumetric  equation 

([/f<°>+£  ffW]  'f41) )  (r) + jv)  (r)  =  -(A'^pf^r)  for  r  €  fi  (C.15) 
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for  the  unknown  fields  pi1*  and  pv,  with  the  r.h.s.  expressed  in  terms  of  the 
previously  determined  surface  problem  solution  pin*. 

The  solution  procedure  amounts  then  to  solving  the  surface  equation 
(C.6),  the  volumetric  equation  (C.15),  and  forming  the  final  solution  ac¬ 
cording  to  Eq.{C.13),  as  described  in  Section  3,2. 

We  stress  the  importance  of  Eq+(C.9)  (he.,  extension  of  Eq,(C,6)  to  the 
interior  of  the  region  SI)  in  the  derivation  of  Eq.(CT5),  Eq,(C.9)  is  necessary 
for  cancellation  of  large  terms  ~  which  would  be  otherwise  present  in 
Eq.(CT5). 


C.2  Solution  of  the  discretized  problem 

In  the  discretized  problem,  the  counterpart  of  the  solution  procedure  de¬ 
scribed  above  is  as  follows: 

In  analogy  to  the  integral  operator  decomposition  (CTO),  we  represent 
the  stiffness  matrix  of  Eq.  (3.6)  as 

A  =  Ar  +  Am  +  Ad  =  ,  (C.16) 

where  the  elements  of  the  matrices  and  arc  given,  explicitly,  by 


A{°1  -  -  y  d:,r !  J  d2r2  [Vr^Q(rx)]  g( rx  -  r2)  ■  n(r2)^(r2)  ,  (C.17a) 


H  tHJ 


and 


=  j  d3rj  j  d2r2  [Vri0a(rl)J  s(rx  -  r2)  ■  n(r2) 

H  an 

+  J  d3r<t>a(  r)^L_0^{r) 

+  j  d3r1d:V20a(r1)g(r1  -  r2)  |  Pa 


*o  rfra)J 

-  /  y  d3r2  tvr0o(ri)]  sfn  -  r2)  ■  ^-)  0,(r2)  . 

(C.lTb) 


n  n\&n 


We  note  that,  when  p(r)/pn  — *  oc  uniformly  in  the  volume  fl  and  £  is  defined 
by  Eq.(CT2),  the  matrices  and  remain  finite. 

It  is  clear  from  Eq.(CT7b)  that  the  matrix  elements  may  be  nonzero 
only  for  basis  functions  <p^  overlapping  the  boundary  dfl  of  the  object.  We 
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call  those  basis  functions  “surface”  basis  functions,  and  the  remaining  ones 
" volume n  basis  functions,  and,  accordingly,  decompose  the  matrix  (C^  16) 
into  blocks  as 


(AS) 

(A& 


(C.18) 


(0) 

where  wo  used  the  fact  that,  by  construction,  the  elements  Aa^  arc  nonzero 
only  for  “surface"  column  indices  0  (eh  Eq,(C.14)). 

Similarly,  we  decompose  the  solution  vector  x  (Eq,(3,3))  into  “surface" 
and  “volume”  blocks  as 


'.T<01  +U'Y 

xv^ 

Xy 

where  we  define  j*s01  as  the  solution  of  the  equation 


(C.  19) 


(C.20) 


here  bs  is  the  surface  block  of  the  r.h.s.  vector  b 


b  = 


bs 

K 


(C.21) 


If  we  now  substitute  Eqs.  (C.18),  (C.19),  and  (C.21)  into  the  original 
equation 

A  x  =  b  ,  (C.22) 

and  use  Eq.(C.2Q),  we  obtain  the  system  of  equations 

(A^  +  £  A™)  s£>  +  xv  =  -Ag  x™  ,  (C.23a) 

(4£  +  €  A$)  x™  +  Ai«  xv  =  -AilJ  xP  +  r1  (K  -  AP  x<°>)  ,  (C.23b) 

where  the  r.h.s. s  arc  expressed  in  terms  of  the  solution  of  Eq.(C,2G). 

As  in  the  derivation  of  Eq.(C.15)  in  the  previous  Section,  we  now  have 
to  show  that  Eq.(C.20)  implies 

4°M0)  =  *>v;  cc.24) 


and  both  these  equations  are  just  the  discrete  counterparts  of  Eqs.  (C.6) 
and  (C.9),  obtained  by  expanding  the  field  into  the  basis  functions  <pa 
(Eq.(3.3))  and  projecting  on  the  same  set  of  basis  functions. 
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It  follows  then  from  Eq.(C.24)  that  Eqs.  (C.23)  reduce  to 


4(1)1 

^VV  _ 


\Awxm 

Ass 

AW  JO) 

/lys  Xs 


=  h 


(C.25) 


We  note  that  the  solution  depends  on  the  original  r.h.s.  6  only  through  its 
surface  block  6S  appearing  in  Eq,(C.20). 

Thus,  finally,  the  solution  procedure  for  the  discretized  equations  is 
the  exact  counterpart  of  the  procedure  described  in  Section  C.l  following 

Eq.(C.15): 

1.  solve  the  surface  problem  (C.20); 

2.  by  using  the  solution  xi°\  compute  the  r.h.s.  of  Eq.(C.25), 

3.  solve  the  volumetric  problem  (C.25)  for  xil]  and  xv;  and 

4.  construct  the  full  solution  according  to  Eq.(C.19). 


C.3  Implementation  in  the  code 

Wc  have  implemented  the  above  two-stage  (surface  +  volume)  solution  pro¬ 
cedure  in  our  code;  some  representative  results  were  described  in  Section  3.4. 
This  approach  enabled  us  to  obtain  accurate  solutions  for  practical  problems 
with  large  p/p0;  when  using  the  original  formulation,  wc  were  not  able  to 
generate  a  meaningful  solution  due  to  instabilities  caused  by  ill-conditioning. 

In  the  present  implementation  of  the  algorithm  wc  construct  the  matrices 
and  >1^  just  by  modifying  the  input  to  the  construction  of  the  full 
matrix  A  Further,  instead  of  solving  Eq.(C.25),  we  solve  an  equivalent 
equation  involving  the  original  matrix  ,4  and  a  “right  preconditioner”  - 
an  additional  diagonal  matrix  (denoted  Q)  which  effectively  rescales  the 
columns  of  the  matrix  A.  Wc  describe  these  element  of  the  implementation 
below: 


Construction  of  the  matrices  and  The  code  uses  as  input 

two  quantities  characterizing  material  properties  of  the  object: 


7*0*) 


*(r) 

kq 


7„(r)  =  1  - 


Po 

p{  r) 


(C.26) 
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In  terms  of  these  functions,  the  matrix  elements  of  Eqs.  (3.6d)  -  (3.6f)  can 
be  written  as 

=  y  d3r^a(r)  [1  -  7p(r}]  ^(r) 

+  k2  j  d3rj  d3r2  0Q(rj)  g{rv  -  r2)  [7«(r2)  +  7p{r2)]  <j>(i( r2)  (C.27) 

+  J d\d3r2  [Vrj0o(r1)]  g(ri  -  r2)  •  [Vrjp(r2)]  00{r2)  . 

It  can  be  now  easily  checked  that  the  matrices  and  of  Eqs,  (C,17) 
are  obtained  from  the  general  equation  (C.27)  with  the  functions  7*  and  7^ 
replaced,  respectively,  by 

7i0)(r)  =  -l,  7f(r)  =  l  (C.28) 

(corresponding  to  /c(r)  =  0  and  p(r)  =  00)  and  by 

7'1)(r)  =  j[l  +  7K{r)]  -1  =  ^-1,  {C.29a) 

<  £*o 

7^(r)  =  1- |[1- 7p(r)]  =  1  -  ■  (C.29b) 


Construction  of  the  rescaled  volumetric  equation* 

matrix  appearing  in  Eq.(C.25)  can  be  represented  as 


We  note  that  the 


AW 

/tsv 

A{1) 

^VV 


L^vs 


€^J  [0 


0 

rlJ 


=  AQ 


(C.30) 


(cf,  Eq.(C.  18)).  Therefore,  instead  of  solving  Eq.(C.25)  we  solve  the  equation 


AQ 


(0.31) 


L^v  J 

involving  the  original  matrix  A  and  the  diagonal  matrix  Q  playing  the  role 
of  a  “right  preconditioncr” . 


D  A  fast  solution  method  for  integral  equations  in 
acoustics 

In  view  of  the  complexity  of  the  anatomical  models  of  the  human  head  and 
the  resulting  problem  sizes  (up  to  millions  of  unknowns),  compression  of  the 
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stiffness  matrix  and  a  fast  solution  method  are  necessities.  In  Phase  I  we 
implemented  an  extension  of  the  FFT-hased  (AIM)  matrix  compression  [2] 
to  acoustics, 

Wc  first  present  the  formulation  of  the  method  in  its  application  to 
acousticsf  and  then  describe  its  implementation  in  our  code. 


D.l  General  formulation 


The  AIM  approach  simply  amounts  to  replacing,  in  computation  of  elements 
of  the  stiffness  matrix,4  the  MoM  basis  functions  (such  as  the  functions  <f>a  in 
Eqs.  (3.3)  and  (3.6))  by  sets  of  “approximately  equivalent”  auxiliary  point 
sources  located  at  nodes  of  a  regular  Cartesian  grid.  The  “equivalence” 
is  understood  in  the  sense  that  the  two  source  representations  generate  the 
same  (to  the  desired  tolerance)  far  field.  Depending  on  the  equation  type,  wc 
apply  the  equivalent  source  expansion  to  the  scalar  basis  functions  mid  their 
gradients  (Eq.(3.6)),  to  individual  components  the  vector  basis  functions, 
say  to  their  divergences,  derivatives  in  the  tensors,  etc. 

To  simplify  the  notation,  wc  denote  in  the  following  by  \p  any  of  those 
basis  functions  or  their  components, 


V>o(r) 


'*.( r), 

V<t>a(r), 

V  ■  tl>a{ r) , 


(D.l) 


In  the  case  of  the  scalar  basis  functions  the  index  a  refers  to  tetrahedra  (for 
constant  basis  functions)  or  vertices  (for  linear  basis  functions),  etc. 

The  compression  procedure  is  as  follows: 

1.  We  construct  a  regular  Cartesian  grid  Q  of  node  spacing  ht  covering 
the  scattcrcr,  and  select  a  multipole  expansion  order  M  =  0, 1,2, — 

2.  Then,  for  each  basis  function  <pa  define  a  cube-shaped  set  Ca  (which 
we  call  the  “expansion  cube”)  consisting  of  (M  +  l)3  Cartesian  grid 
nodes  u  £  S,  and  covering5  in  the  optimal  way  (according  to  predefined 

4 More  precisely,  the  replacement  is  done  in  the  terms  involving  the  Green  function. 

HVe  have  to  require  here  that  Mh  is  larger  than  than  the  diameter  of  the  basis  function 
support. 
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geometrical  criteria)  the  basis  function  support.  For  example,  we  may 
require  that  the  center  cQ  of  the  cube  CQ  is  located  closest  to  the 
centroid  of  the  basis  function  support  (the  center  cQ  would  coincide 
with  one  of  grid  nodes  for  M  even,  and  with  a  grid  cell  center  for  M 
odd). 

3.  We  approximate  the  original  basis  function  by  a  set  of  equivalent 
point  sources  at  the  nodes  of  the  expansion  cube  Ca, 

*  &(r)  =  £  Vau  <S3{r  -  u)  ,  (D.2) 

ueCa 

where  the  set  of  (M  +  I  )3  coefficients  V  (strengths  of  the  sources)  is  de¬ 
termined  so  that  both  functions  have  the  same  the  multipole  moments 
up  to  the  order  M,  i.eM 

j  d 3ppm  [<pa(cQ  +  p) -£a{ca  +  p)] 

=  J  d3p  p?T  pyv  pTz  +  p)  -  0a(cQ  +  p)) 

=  0  for  0  <  mx,  my,  mz  <  M  . 

Here  and  below  we  use  the  multi-index  notation  with  m  =  (mx  my  m2), 
and  the  multipole  moments  arc  evaluated  relative  to  the  center  cQ  of 
the  expansion  cube* 

In  order  to  find  the  set  of  (M  -f  l)3  coefficients  V 

(a)  We  evaluate  the  set  of  (M  +  l)3  multipole  moments  up  to  the 
order  M, 

(<P<>) m  (cn )  =  f  d3p  pm  <pa (ca  +  p)  for  0  <mx,mv,mt<  M 

(DA) 

of  the  basis  functions  with  respect  to  the  expansion  cube  center 

cc*- 

(b)  We  solve  the  set  of  (M  +  l)3  equations 

Y  (u  -  ca)“  vnil  =  (ipa)m( c0)  for  0  <  mx, my,  mt  <  M  . 

uecn 

(D.5) 
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Eq,(D,5)  is  the  Vandermonde  system,  and  its  well-known  closed- 
form  solution  is  given,  e.g.,  in  Ref. [2]  and  in  Appendix  B.2  of  our 
Status  Report  2. 


4.  In  the  near-field  range  wc  compute  the  matrix  elements  according  to 
the  original  MoM  expression  (such  as  (3*6)),  while  in  the  far-ficid  range 
wc  utilize  the  expansion  (D,2)  for  the  basis  functions.  This  yields  the 
representation  of  the  matrix  dements 


Near  .  Far 
a<*0-a<*0  +aaP 


where 


Near  _ 

- 


aad  ~  ^or  dist(Ca  —  Cp)  <  dh  and 


Ka  ~  Q^rl 

laa/3l 


otherwise 


and 


(D.G) 

>  e  , 
(D.7) 

CD-8) 


£0=  E  ^us(u-v)v/Jv. 

u 

u^v 

In  Eq.(D.7)  dist(«  *  * )  denotes  the  distance  between  the  expansion  cubes, 


dist(Cft  —  Cp)  =  min  ||u  —  v||  (D*9) 

u  eCo.veCo 


and  the  distance  [|  ■  ||  (in  fact,  the  ||  ■  norm)  is  defined  by 

[lull  =  maxfluj,  |u„|,  |»s|}  .  (D.10) 

The  integer  parameter  d  in  Eq.(D.7)  is  the  near- field  range ,  and 
the  parameter  e  >  0  is  the  near- field  tolerance .  Eq.(D.7)  states 
that  in  the  near-field  range  the  matrix  element  aaii  is  replaced  by  the 
difference  of  its  exact  value  and  the  far-field  approximation  «3f. 
less  the  approximation  is  accurate  to  the  tolerance  £,  The  subtraction 
of  the  far-field  term  is  necessary,  since,  for  arbitrary  locations  of  the 
basis  function  supports,  there  is  no  translationally  invariant  way  of 
expressing  the  geometrical  ncar-ficld  condition  in  Eq.(D.7)  6  -  and 
the  factorized,  translationally  invariant  form  of  Eq.(D.S)  is  the  crucial 

6 We  note  that,  on  the  other  hand,  the  restriction  u  /  v  in  Eq.fD.8)  is  equivalent  to 
the  translationally  invariant  condition  j?(0)  =  0. 
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element  of  the  representation.  In  fact,  the  far-field  part  (D.8)  of  the 
matrix  is  a  product  of  the  Toeplitz  matrix  G  of  elements  g(u  -  v}T  and 
the  sparse  matrix  V  and  its  transpose;  clearly,  the  Toeplitz  property 
allows  computation  of  the  matrix-vector  product  by  means  of  Fast 
Fourier  Transform  (FFT)*  Since,  by  construction,  the  near-field  part. 
(D.7)  of  the  matrix  is  sparse,  Eq.(D.6)  provides  a  compressed  repre¬ 
sentation  of  the  original  matrix. 

For  the  given  geometry  and  its  MoM  discretization,  the  error  due  to  the 
matrix  compression  (Section  D.4)  and  the  computational  cost  (Section  D.6) 
are  controlled  by  the  parameters  h ,  d,  and  M ,  be.,  the  Cartesian  grid  spac¬ 
ing,  the  near-field  range,  and  the  multipole  expansion  order.  We  discuss 
later  on  the  optimal  choice  of  these  parameters,  ensuring  the  desired  accu¬ 
racy  and  minimizing  the  computational  cost. 

D.2  Application  to  specific  integral  operators  in  acoustics 

Wc  present  here  applications  of  the  AIM  compression  to  some  specific  inte¬ 
gral  operators  and  their  discretized  representations* 

Constant  basis  functions.  In  this  case  we  apply  AIM  compression  to 
the  monopole  and  dipole  matrices  of  Eqs,  (B.l).  In  the  first  of  these  wc  use 
the  expansion 

<Mr)  *  £,<r)  =  Y,  *3<r  -  u)  (D  U) 

for  the  basis  functions  constant  on  tetrahedra  a,  where  the  expansion  coef¬ 
ficients  Vu  arc  evaluated  as  described  in  Appendices  B.l  (evaluation  of  the 
multipole  moments)  and  B.2  (solution  of  the  Vandermonde  system  (D<5)) 
of  our  Status  Report  2, 

Wc  obtain  in  this  way  the  far-field  approximation 

(“"“I  •  (D.12) 

u,v  \K0  p0) 

Similarly,  in  the  dipole  term  wc  approximate  the  gradient  of  the  basis 
function  and  its  product  with  the  material-dependent  function  A 

V0Q(r}*  E  v°u53(r-u)  (D.13a) 

uec0 
and 
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(D.13b) 


V0,i(r)  A(r)  ~  Yi  v«u<S3(r  -  u) 

uecQ 


Wo  find  in  this  way 


4^EV^(u-  v)-V°v  >  (D-14) 

UfV 

where  the  vector- valued  expansion  coefficients  VI}  and  VDA  arc  given  by 


V‘J„  a  y  n,  V, 

aU  ^  fa  fa  U 

(D.15a) 

v™=  E  */.  a/„v/„u! 

(D.  15b) 

/ 


here  V ju  are  expansion  coefficients  for  a  constant  basis  function  supported 
on  the  triangular  face  /,  the  sums  are  taken  over  all  faces  of  the  consid* 
ered  tetrahedra,  We  note  that  the  expansion  coefficients  VDA  depend  on 
the  discontinuities  A  of  the  material  parameter  p^jp  on  the  individual  faces. 
In  this  way,  Eqs,  (DT2)  and  (D.I4)  provide  a  compressed  matrix  represen¬ 
tation  on  terms  of  scalar  expansion  coefficients  VM  and  two  sets  of  vector 
expansion  coefficients,  VD  and  VDA  (hence,  the  total  of  7  components). 
Explicit  expressions  for  the  vectorial  expansion  coefficients  can  lie  obtained 
by  computing  the  inultipolc  moments  of  gradients  of  the  basis  functions  and 
solving  the  resulting  Vandermonde  system.  Such  results  arc  given  in  Rcfi  [2]; 
we  also  provided  specific  expression  for  the  acoustics  problems  in  our  Status 
Report  2. 

We  also  note  that  the  expansion  coefficients  V^u  of  Eq,(DT5)  represent 
the  gmdient  V0a  of  the  basis  function  Therefore,  the  zero-th  order 
(monopolc)  of  the  far-ficld-equivalent  basis  function  (D.13a)  must  vanish 
exactly;  this  requires  a  high  numerical  precision  in  evaluating  the  coefficients 
VD  and  the  pertinent  far-field  expansions.  A  possible  alternative  way  of 
parameterizing  the  F FT- com  pressed  matrix  could  be  to  use  dipole  rather 
than  monopolc  point  sources,  which  would  then  automatically  guarantee 
vanishing  of  the  monopole  component  of  the  far- field.  On  the  other  hand,  the 
expansion  (DT3b)  will,  generally,  yield  a  nonvanishing  monopolc  moment, 
and  will  be  easier  to  handle  numerically 

Linear  basis  functions.  The  structure  of  the  basis  function  and  matrix 
element  expansions  is  here  the  same  as  in  Eqs*  (DTI)  to  (D.15b)*  Explicit 
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expressions  for  the  expansion  coefficients  can  be  obtained  by  computing 
the  multipole  moments  of  gradients  of  the  basis  functions  as  described  in 
Appendix  B.1,3  of  our  Status  Report  2, 


D.3  Geometry  and  discretization 

Since  we  are  concerned  with  acoustic  wave  propagation  in  biological  media, 
we  concentrate  on  volumetric  problems,  assuming  an  approximately  uniform 
distribution  of  the  material  in  the  considered  volume*  For  the  purpose  of 
algorithm  complexity  estimates,  we  introduce  the  following  quantities  and 
notation: 

1.  Wc  consider  an  object  {e.g.,  a  cube)  of  linear  sizes  L  x  L  x  L  uni¬ 
formly  discretized  with  the  spatial  resolution  a *  The  number  of 
MoM  unknowns  is  then 


(D.16) 


with  a  coefficient  cv  dependent  on  the  discretization*  E.g*,  for  the 
case  of  partition  into  regular  cubes,  with  each  of  them  partitioned 
into  a  minimum  number  of  tetrahedra,  wc  have  cv  =  5.  Wc  stress 
that,  in  this  context,  the  parameter  a  controls  merely  the  number  of 
unknowns  per  unit  volume  of  the  scatterer.  Wc  only  assume  that  the 
sizes  of  basis  function  supports  are  not  larger  than  a. 

2*  We  introduce  a  regular  Cartesian  grid  Q  of  K  x  A"  x  A"  nodes  with  the 
cell  size,  or  grid  spacing  /i,  hence  A'  =  L/h.  Therefore,  the  number 
of  grid  nodes  is 


(D.  1 7) 


Wc  allow  the  grid  spacing  h  to  be  large  relative  to  the  geometry  spatial 
resolution  6;  however,  we  will  find  that  the  optimal  grid  spacing  is 
comparable  to  a. 

3.  With  the  partition  specified  above,  the  number  of  unknowns  nc  per 
grid  cell  is  then,  on  the  average, 


(D.  18) 
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D.4  Error  estimates 


Since  we  arc  mostly  concerned  here  with  subwaveiengtb  problems,  we  can 
determine  approximate  error  bounds  by  neglecting  the  oscillatory  behavior 
of  the  Green  function  (3*2)  altogether,  i.e.,  considering  the  Laplace  equation 
Green  function 

(D19) 

The  estimates  wc  obtain  remain  valid  as  long  as  the  basis  function  supports 
and  the  Cartesian  grid  spacing  h  are  small  compared  to  the  wavelength.  Wc 
use  the  assumptions  about  the  geometry  and  its  discretization  as  specified 
in  Section  D.3* 

The  error  in  the  matrix  elements  due  to  the  approximation  (D*2)  is 
controlled  by  the  (integer)  near- field  range  d  and  the  multipole  expansion 
order  M.  According  to  the  definition  of  the  range  d  in  Eqs.  (D.7)  and  (D.9), 
the  minimum  far-field  distance  between  the  centers  of  the  expansion  cubes 
is 

M)  =  (rf  +  M)  h  (D.20) 

(wc  recall  that  the  cube  size  is  M  h).  A  qualitative  estimate  of  the  error  in 
the  matrix  elements  due  to  the  multipole  approximation  (D.2)  is 


(  1  \M+l 

\d  +  M  ) 


(D.21) 


However,  since  the  coefficient  in  this  estimate  is  not  easy  to  determine  ana¬ 
lytically,  we  use  numerical  data. 

In  Fig.  19  we  show  the  computed  relative  errors  in  matrix  elements 
for  scalar  piecewise  constant  functions  0Q  (Appendix  B.i),  computed  for 
expansion  orders  M  =  2,3,4,  and  for  grid  spacing  h  equal  to  the  basis 
function  support  size  a.  The  errors  are  plotted  as  functions  of  the  distance 
R  between  the  expansion  cubes  (measured  in  the  grid  spacing  units  h ),  and 
the  expansion  order  M,  for  typical  values  of  the  grid  spacing  relative  to 
the  basis  function  support  sizes.  The  plot  is  based  on  a  statistical  sample  of 
matrix  elements  for  various  distances  and  orientations  of  the  basis  functions1 
supports.  The  errors  depend  somewhat  on  the  ratio  of  the  grid  spacing 
h  to  the  support  size  a,  but  remain  practically  unchanged  in  the  range 
|  a  <  h  <  2  a. 
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AIM  errors  -  scalar  basis  functions 


Figure  19:  Relative  errors  of  the  AIM  approximation  to  matrix  elements  for 
scalar  basis  functions,  for  expansion  orders  M  ~  2, 3, 4,  plotted  as  functions 
of  the  distance  R  between  the  celt  centers.  The  lines  indicate  approximate 
“envelopes"  of  maximum  errors. 

The  Table  4  summarizes  results  of  the  error  analysis  in  terms  of  the 
required  near-field  range  as  a  function  of  the  expansion  orders  M  and  the 
desired  accuracy  in  the  matrix  elements,  from  LO  %  to  0.2  %♦  The  parameter 
d  in  the  Table  (defining  the  near- field  range  in  AIM  through  Eq.(D>20))  were 
determined  on  the  basis  of  numerically  computed  errors,  such  as  shown  in 
Fig.  19. 


Table  4:  AIM  near  field  range  d  as  a  function  of  the  expansion  order  M  and 
error  tolerance 


AIM:  d 

M 

scalar  basis  functions 

1.0% 

0.5% 

0.2% 

2 

2 

3 

4 

3 

1 

1 

2 

4 

0 

1 

1 
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D*5  Near-  and  far- field  computation  in  AIM 

We  present  here  the  main  elements  of  the  implementation  of  the  AIM  algo¬ 
rithm  for  matrix-vector  multiplication;  they  will  be  used  in  the  following  to 
obtain  the  computational  complexity  of  the  algorithm* 

In  the  compressed  matrix  representation  (D*6)  the  near-field  part  of  the 
matrix  (D,7)  is  the  conventional  MoM  matrix  restricted  to  the  near-field 
range  and  stored  in  a  sparse  form.  Matrix- vector  multiplication  involving 
this  matrix  is  implemented  in  straightforward  way  as  a  sparse-mat rix- vector 
product. 

Eq.(D.8)  implies  that  the  far-ficld  part  of  the  mat  rix- vector  product  is 

N 

?/oar  =  5I  S  Vau3(u  "  v)  Vpvxg  (D.22) 

0=1  utv£Q 

where  the  sums  arc  taken  over  all  Cartesian  grid  nodes  (cell  centers)  u  and 
v. 

Equivalently,  the  matrix-vector  multiplication  (D.22)  can  be  executed  as 
the  sequence  of  three  steps 

N 


(A): 

£V  Vl*vX0 

0=1 

for  all  v  6  Q  , 

(D.23a) 

(B): 

Vu  =  53  -  v) 

for  all  u  €  G  r 

(D.23b) 

(C): 

vl"=  X]  Vaufiu 

U€Q 

for  all  o  =  1, 2, .  * .  ,  N  . 

(D.23c) 

The  convolution  over  the  Cartesian  grid  nodes  in  Eq.(D/23b)  is  imple¬ 
mented  by  means  of  an  FFT:  We  compute  the  FFTs  of  the  Green  function 
expansion  coefficients, 


(BO):  ff(Q)= 

for  Q  €  0  , 

(D.24a) 

vzG 

of  the  input.  Cartesian  vector  x, 

(Bl)  :  ^  FQvxv 

for  Q  €  G  , 

(D.24b) 

veG 

(where  the  matrix  T  represents  the  discrete  Fourier  transform,  and  Q  is  the 
conjugate  Cartesian  grid),  we  evaluate  the  product  of  the  Fourier  transforms, 
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{ D .  24c ) 


(B2)  :  yQ  =  g{ Q)  xQ  for  Q  e  G  , 

and  take  the  inverse  Fourier  transform  of  the  product, 

(B3)  :  ?u  -  ^  -^uQ  Vq  for  u  ^  &  ■ 

QtG 


(D.24d) 


Of  course,  the  FFT  of  the  Green  function  in  step  (BO)  lias  to  be  computed 
and  stored  only  once. 

D.6  Computational  complexity  in  AIM 

The  computational  complexity  estimates  we  give  here  are  obtained  for  a  vol¬ 
umetric  problem,  with  the  simple  geometry  and  Cartesian  grid  as  specified 
before  in  Section  D<3. 

Near- field  computation.  According  to  Eq*(D.20),  the  AIM  near-field 
volume  of  a  given  basis  function  is  a  cube  of  side  length  2(d  +  M)  h}  con¬ 
taining  on  the  average  nc  (Eq.(D.lS))  unknowns.  It  follows  that,  the  total 
number  of  near-held  matrix  dements  is 


=  5  N  [2  (d  +  M )]3  nc  =  4  (rf  +  M)3  —  , 

Is  G 


(D.25) 


with  the  factor  \  due  to  the  symmetry  of  the  matrix.  Consequently,  the 
near-field  computation  requires 

N2 


Nop  [near]  A]M  =  8  {d  +  M) 


N, 


(D.26) 


operations  (complex  multiplications  and  additions). 


Far-field  computation.  Conversion  of  the  set  of  MoM  unknowns  to  the 
Cartesian  vector  in  step  (A)  (Eq.(D,23a))  involves  (M  + 1)3  grid  points  in  an 
expansion  cube  for  each  unknown;  the  same  applies  to  the  inverse  operation 
(C)  of  Eq.(D.23c).  The  number  of  complex  multiplication  and  additions  is 
thus 


Nop  [step  A]aim  -  Nop  [step  CjA1M  =2  (M  +  1)‘*  N  . 

In  the  FFT-bascd  implementation  of  step  (B),  he.,  steps  {Bl} 
Eqs.  (D/24),  the  numbers  of  operations  arc 

Nop[stcp  B1]aim  =  Nop  [step  B3]AIM  =  \  SNg  log2(8  Ng) 

=  20  Ng  log  2  (8 Ng)  , 


(D.27) 
(B3)  of 

(D.28) 


47 


where  the  factor  8  comes  from  the  zero  padding  of  the  Cartesian  vectors, 
and  where  wc  assumed  the  FFT  complexity  as  for  the  radix-2  algorithm, 
i.e.,  |  n  log2n  for  a  vector  of  length  n.  Step  (BO),  which  has  to  be  carried 
out  only  once,  belongs  to  the  compressed  matrix  set-up,  and  we  don't  count 
it  as  a  part  of  matrix- vector  multiplication. 

Finally,  step  (B2)  of  Eq+(D,24c)  -  multiplication  of  Cartesian  vectors  - 
requires 

Nopfstep  B2]A1M  =  WNg  (D.29) 

operations.  The  total  numbers  of  operations  in  the  far-ficld  computation 
(steps  (A),  (B),  and  (C))  is  thus 

Nop[far]AIM  =  4  (M  +  l)3  N  +  8  [2  +  5  log2(8 JVC)]  Ng  .  (D.30) 

Total  computational  cost*  By  collecting  contributions  of  Eqs,  (D.26) 
and  Eqs>  (D,30)  wc  find  the  total  computational  complexity  of  matrix- vector 
multiplication  in  AIM, 


kt2 

Nop  [A  *]AW  =8  (d+  A/)3  —  +  4  {M  +  l)3  N  +  8  [2  +  5  log  2(8Ne)]  Ng 


Thus,  if  we  define 

i 

III 

Si55 

(D.31) 

(D.32) 

the  computational  cost 

per  unknown  is 

—  Nop[AxJAIM  = 

^ +a2(M)  +  [a3  +  a4  ln{£JV)]g 

(D.33) 

with 

ax{d,M)  =  8{d  +  Mf  , 

(D.34a) 

o2(M)  =  4(M+1)3  , 

(D.34b) 

a3  =  136  , 

(D.34c) 
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a4  =  —  -  577  . 

4  In  2 

(D.34d) 

We  list  in  Table  5  the  values  of  the  coefficients  for  several  choices  of  the 
near-field  ranges  and  expansion  orders,  based  on  the  error  estimates  sum¬ 
marized  in  Table  4.  Wc  have  selected  the  parameters  { d ,  M)  corresponding 
to  the  expected  1%  error  level,  as  given  in  Table  4.  In  the  ease  of  AIM  we 
chose  the  parameters  for  the  less  favorable  case  of  vector  basis  functions. 
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Table  5:  Coefficients  in  matrix- vector  multiplication  algorithm  complexity 
in  AIM 


AIM 

d 

M 

a  M) 

a2{M) 

a4 

10 

2 

13,824 

108 

136 

57.7 

3 

3 

5,532 

256 

136 

57.7 

D.7  Storage  in  AIM 

Storage  (in  complex  numbers)  required  for  the  compressed  matrix  represen¬ 
tation  in  the  AIM  algorithm  is  obtained  in  close  analogy  to  the  complexity 
estimates: 

A j2 

Mem  [A*™1  +  /tFar]AIM  =  4  (d  +  Mf  —  +  {M  +  l)3  N  +  8Ng  ,  (D.35) 

NG 

i.e*,  in  terms  of  the  ratio  £  (Eq.(D.32))  the  total  storage  per  unknown  is 

^  Mem[j4]AIM  =  ^Ml+a2(M)  +  a.^  (D.36) 


with 


al(d,M)  =  4(d+M)3  ,  (D,37a) 

a2{M)  =  (M  +  l)3  ,  (D.37H) 

a3  —  8  .  (D*37c) 

The  values  of  the  coefficients  2-  for  typical  parameters  d  and  M  arc  given 
in  the  Table  6. 

Table  fi:  Coefficients  in  relative  matrix- vector  matrix  storage  in  AIM 


AIM 

d 

M 

3,  KM) 

S  2  (AO 

aA 

10 

2 

6,912 

27 

8 

3 

3 

2,766 

64 

8 
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D.8  Optimization  of  parameters  in  AIM 

In  Figs.  20  and  21  we  plot  the  computational  costs  and  memory  require¬ 
ments  (per  unknown)  for  the  AIM  algorithm  in  a  rather  wide  range  of  the 
parameter  £  —  NgfN>  covering  typical  values  encountered  both  in  volume 
and  in  surface  problems. 


Figure  20:  The  matrix- vector  multiplication  cost-per-unknown  in  the  AIM 
algorithm  for  N  —  10°,  plotted  as  a  function  of  £  =  Ng/Nf  in  the  doubly* 
logarithmic  scale. 


Figure  21:  The  matrix  storage  per  unknown  in  the  AIM  algorithm,  plotted 
as  a  function  of  £  =  Ng/N ,  in  the  doubly- logarithmic  scale. 

Figs.  20  and  21  show  that  the  minimum  cost  of  the  matrix- vector  mul¬ 
tiplication  (which  we  take  as  the  decisive  factor  in  the  optimization)  is 
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achieved  at  £  =  Ng/N  oc  2.  According  to  the  discussion  of  Section  D+3, 
the  grid  spacing  corresponding  to  this  value  of  £  is 


h  o*  a  (cy  £)  ^  0.46  a  , 


(D.38) 


where  a  is  the  spatial  resolution  of  the  MoM  discretization.  Clearly,  this 
value  should  not  be  taken  too  literally,  since  the  actual  number  of  MoM 
unknowns  may  strongly  depend  on  the  geometry;  nevertheless,  Eq.(D.38) 
indicates  that  the  optimal  grid  spacing  is  close  to  (and  likely  smaller 
than)  the  resolution  of  the  geometry  discretization. 

Fig.  21  indicates  that  the  minimum  storage  is  attained  for  a  higher  value 
of  f  ^  20,  corresponding  to  the  grid  spacing  smaller  by  another  factor 
IQ-1/3  ^0.46.  The  estimated  storage  at  £  =  2  is  exceeds  its  minimum  value 
by  about  a  factor  4.  Therefore,  if  memory  reduction  is  of  higher  priority,  the 
grid  spacing  may  be  decreased  accordingly  to  reach  the  overall  optimum. 

We  also  see  from  Fig.  20  that  the  quadrupole  expansion  (order  M  —  2) 
yields  the  cost  almost  a  factor  of  2  lower  than  the  next  order  (M  =  3) 
expansion.  At  the  same  time,  Fig.  19  shows  that,  for  acoustics,  the  near- 
field  range  of  3  to  4  grid  spacings  is  sufficient  to  reduce  the  far-fteld 
error  to  the  1%  level.  Further  empirical  analysis  has  shown  that  the  optimal 
compression  parameters ,  ensuring  errors  not  larger  than  1  %,  are 


*  -  2a  >  -3  h  ,  M  =  2 


(D.39) 


(cf,  Eq.(D.20)).  This  conclusion  applies  to  problems  with  moderately  to 
highly  subwavclcngth  discretization,  i.c.,  with  tetrahedron  size 


A 


(D.40) 


i.e.,  in  practice,  to  all  problems  we  arc  considering  here;  it  has  been  con¬ 
firmed  by  many  practical  computations  for  large  problems  (up  to  several 
million  unknowns)  which  wc  report  on  in  Section  3.4. 
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